#include "SignalGeneration.hpp"

#include <math.h>
#include <time.h>

#include <string.h>
#include <stdlib.h>

#include <stdexcept>
#include <string>

using namespace std;

static const double PI = atan(1.0) * 4;



//
// CSignalGenerationInfo
//

CSignalGenerationInfo::CSignalGenerationInfo(double TotalObservationTime, int SampleCount)
	: fTotalObservationTime(TotalObservationTime), fSampleCount(SampleCount), fNoiseAmplitude(0), fPhase(0), fOscillatorCount(0)
{
}

void CSignalGenerationInfo::AddOscillator(double Frequency, double Amplitude)
{
	Frequencies.push_back(Frequency);
	Amplitudes.push_back(Amplitude);
	fOscillatorCount += 1;
}

int CSignalGenerationInfo::OscillatorCount()
{
	return fOscillatorCount;
}

int CSignalGenerationInfo::SampleCount()
{
	return fSampleCount;
}

double CSignalGenerationInfo::NoiseAmplitude()
{
	return fNoiseAmplitude;
}

void CSignalGenerationInfo::NoiseAmplitude(double value)
{
	fNoiseAmplitude = value;
}

double CSignalGenerationInfo::TotalObservationTime()
{
	return fTotalObservationTime;
}

double CSignalGenerationInfo::GetFrequency(int Oscillator)
{
	if (Oscillator > fOscillatorCount)
	{
		throw invalid_argument("Oscillator is out of range");
	}
	return Frequencies[Oscillator];
}

double CSignalGenerationInfo::GetAmplitude(int Oscillator)
{
	if (Oscillator > fOscillatorCount)
	{
		throw invalid_argument("Oscillator is out of range");
	}
	return Amplitudes[Oscillator];
}

double CSignalGenerationInfo::Phase()
{
	return fPhase;
}

void CSignalGenerationInfo::Phase(double value)
{
	fPhase = value;
}



//
// CSamplingGenerationInfo
//

CSamplingGenerationInfo::CSamplingGenerationInfo(double Burst0DeltaT, int PulsesPerBurst)
	: fBurst0DeltaT(Burst0DeltaT), fPulsesPerBurst(PulsesPerBurst), fNumberOfBursts(0)
{
}

void CSamplingGenerationInfo::AddBurst(double StaggerFactor, double BurstPauseTime)
{
	StaggerFactors.push_back(StaggerFactor);
	BurstPauseTimes.push_back(BurstPauseTime);
	fNumberOfBursts += 1;
}

int CSamplingGenerationInfo::NumberOfBursts()
{
	return fNumberOfBursts;
}

double CSamplingGenerationInfo::GetStaggerFactor(int Burst)
{
	if (Burst > fNumberOfBursts)
	{
		throw invalid_argument("Burst is out of range");
	}
	return StaggerFactors[Burst];
}

double CSamplingGenerationInfo::GetBurstPauseTime(int Burst)
{
	if (Burst > fNumberOfBursts)
	{
		throw invalid_argument("Burst is out of range");
	}
	return BurstPauseTimes[Burst];
}

int CSamplingGenerationInfo::PulsesPerBurst()
{
	return fPulsesPerBurst;
}

double CSamplingGenerationInfo::Burst0DeltaT()
{
	return fBurst0DeltaT;
}



//
// CSamplingInfo
//

CSamplingInfo::CSamplingInfo(vector<double> &TimeVector, vector<double> &SampleFrequencies, double TotalObservationTime, vector<double> &DeltaT, int TotalSampleCount)
	: fTimeVector(TimeVector),
	SampleFrequencies(SampleFrequencies),
	fTotalObservationTime(TotalObservationTime),
	DeltaT(DeltaT),
	fTotalSampleCount(TotalSampleCount)
{
}

double CSamplingInfo::TotalObservationTime()
{
	return fTotalObservationTime;
}

double CSamplingInfo::GetSampleFrequency(int Burst)
{
	if (Burst > (int)SampleFrequencies.size())
	{
		throw invalid_argument("Burst is out of range");
	}
	return SampleFrequencies[Burst];
}

vector<double> CSamplingInfo::TimeVector()
{
	return fTimeVector;
}

double CSamplingInfo::TimeVector(int Index)
{
	return fTimeVector[Index];
}

int CSamplingInfo::TotalSampleCount()
{
	return fTotalSampleCount;
}

double CSamplingInfo::GetDeltaT(int Burst)
{
	if (Burst > (int)DeltaT.size())
	{
		throw invalid_argument("Burst is out of range");
	}
	return DeltaT[Burst];
}



// 
// CSignal
//

CSignal::CSignal(vector<complex<double> > &Signal, vector<double> &TimeVector, int SampleCount, double TotalObservationTime)
	: fSignal(Signal), fTimeVector(TimeVector), fSampleCount(SampleCount), fTotalObservationTime(TotalObservationTime)
{
}

vector<complex<double> > CSignal::Signal()
{
	return fSignal;
}

vector<double> CSignal::TimeVector()
{
	return fTimeVector;
}

int CSignal::SampleCount()
{
	return fSampleCount;
}

double CSignal::TotalObservationTime()
{
	return fTotalObservationTime;
}

double CSignal::GetDeltaT()
{
	return fTotalObservationTime / fSampleCount;
}



//
// GenerateSampling
//

CSamplingInfo CSamplingGenerationInfo::GenerateSampling()
{
	double *StaggerRatio = new double[fNumberOfBursts];

	for (int i = 0; i < fNumberOfBursts; i++)
	{
		StaggerRatio[i] = StaggerFactors[i] / StaggerFactors[0];
	}

	vector<double> TimeVector;
	vector<double> SampleFrequencies;
	vector<double> DeltaT;

	double Now = 0;
	for (int i = 0; i < fNumberOfBursts; i++)
	{
		double CurrentBurstDeltaT = StaggerRatio[i] * fBurst0DeltaT;
		DeltaT.push_back(CurrentBurstDeltaT);
		SampleFrequencies.push_back(1.0 / CurrentBurstDeltaT);

		for (int j = 0; j < fPulsesPerBurst; j++)
		{
			TimeVector.push_back(Now);
			Now += CurrentBurstDeltaT;
		}

		Now += -CurrentBurstDeltaT + BurstPauseTimes[i];
	}

	int TotalSampleCount = fNumberOfBursts * fPulsesPerBurst;
	return CSamplingInfo(TimeVector, SampleFrequencies, Now, DeltaT, TotalSampleCount);
}



//
// GenerateSignal
//

CSignal CSignalGenerationInfo::GenerateSignal()
{
	vector<complex<double> > Signal;
	vector<double> TimeVector;
	double CurrentTime = 0;

	// SampleCount - 1 -> Last sample aligns *exactly* at the same time as in the sampled time vector
	double DeltaT = fTotalObservationTime / (fSampleCount - 1);

	for (int i = 0; i < fSampleCount; i++)
	{
		double Re = 0;
		double Im = 0;

		// base signal (no noise)
		for (int j = 0; j < fOscillatorCount; j++)
		{
			double Arg = 2.0 * PI * Frequencies[j] * CurrentTime + fPhase;
			Re += Amplitudes[j] * cos(Arg);
			Im += Amplitudes[j] * sin(Arg);
		}

		// adding noise
		// trigonometric argument of a random noise frequency
		double RandArg = 2.0f * PI * ((double)rand()) / (double)RAND_MAX;
		// the random frequency's amplitude
		double RandAmp = ((double)rand()) / ((double)RAND_MAX / 2.0) - 1.0;
		// add (weighted by noiseamp)
		Re += fNoiseAmplitude * RandAmp * cos(RandArg);
		Im += fNoiseAmplitude * RandAmp * sin(RandArg);
		
		Signal.push_back(complex<double>(Re, Im));
		TimeVector.push_back(CurrentTime);
		CurrentTime += DeltaT;
	}

	return CSignal(Signal, TimeVector, fSampleCount, fTotalObservationTime);
}



//
// GetSignalSamples
//

void CSignal::GetSignalSamples(double *Re, double *Im, CSamplingInfo &SamplingInfo)
{
	if (SamplingInfo.TotalSampleCount() > fSampleCount)
	{
		throw logic_error("The generated signal cannot be sampled because the sampling is too fast.");
	}

	if (SamplingInfo.TotalObservationTime() > fTotalObservationTime)
	{
		throw logic_error("The generated signal cannot be sampled because it is too short time-wise.");
	}

	vector<double> SampledTimeVector = SamplingInfo.TimeVector();

	int SignalIndex = 0;
	for (int i = 0; i < SamplingInfo.TotalSampleCount(); i++)
	{
		double DestinationTime = SampledTimeVector[i];
		
		// loop until the signal time passes the destination time
		while (fTimeVector[SignalIndex] < DestinationTime && SignalIndex < fSampleCount - 1)
		{
			SignalIndex += 1;
		}

		// now the signal pointer 'signal->t[signal_index]' is right after
		// or exactly on top of the destination time.

		// just take the last sample from the input signal. if the original signals' sample rate
		// is sufficient, this is accurate enough.
		Re[i] = (double)fSignal[SignalIndex].real();
		Im[i] = (double)fSignal[SignalIndex].imag();
	}
}


